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ABSTRACT 

We present single-dish 350 fim dust continuum polarimetry as well as HCN and 
HCO + J = 4 — » 3 rotational emission spectra obtained on NGC 1333 IRAS 4. The 
polarimetry indicates a uniform field morphology over a 20" radius from the peak con- 
tinuum flux of IRAS 4A, in agreement with models of magnetically supported cloud 
collapse. The field morphology around IRAS 4B appears to be quite distinct however, 
with indications of depolarization observed towards the peak flux of this source. Inverse 
P-Cygni profiles are observed in the HCN J = 4 — * 3 line spectra towards IRAS 4A, 
providing a clear indication of infall gas motions. Taken together, the evidence gathered 
here appears to support the scenario that IRAS 4A is a cloud core in a critical state of 
support against gravitational collapse. 

Subject headings: ISM: individual (NGC 1333 IRAS 4) — ISM: molecules — polarization 
— submillimeter — dust — stars: formation 

1. Introduction 

The process of star formation is a key phenomenon in astrophysics. It touches upon a diverse 
range of topics including galactic evolution, stellar evolution, planet formation, and astrobiology. 
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Yet d espite the fundamental nature of this process, it remains poorly understood (jCrutcher et al 
2008). A key issue is determining the support mecha nism(s) that gov erns star formation, as it is 
clear stars cannot be forming at a free fall timescale (|Shu et al.lll987l ). Two competing ideas cur- 
rently strive to explain the support mechanism; mag netism plus ambipolar diffusion and turbulence 
plus a weak magnetic field (|McKee Ostrikerll2007l ). At present, a consensus on the nature of the 
support mechanism has not been realized. 

A particularly important quantity for star formation theor y is the mass-to-flux r atio, M/$, 
where <I> is the magnetic flux within a region enclosing a mass M ([Nakamura &: Lill2008l). Magnetic 
supp ort models predict the critical value of this ratio to be eq ual to c$/y/G (Mouschovias &: Spitzer 



19761 ). where G is the gravitational constant and c$ = 0.12 (jTomisaka et al,lll988l ). It follows that 



the value of Mj $ normalized toc$/vG should be close to unity in regions near the point of collapse. 
Turbulent models place no such restrictions on the mass-to-flux ratio but do predict a chaotic field 
morphology assuming that flux-freezing holds. The mass-to-flux ratio can be determined from 
observables and thus may be used as an important indicator as to the nature of the support 
mechanism. 



At a distance of ~ 300 pc0 ( Girart et al. 2006 ). a known site of clustered low an d intermediate 
mass star formation, and possessing young embedded cores at an age of ~ 1 Myr (jHatchell et al 



20051 ). NGC 1333 is an ideal target for studying the onset of clustered star formation. Two such em- 
bedded cores of interest include NGC 1333 IRAS 4A (a 2 ooo = 3h29ml0.42s, (5 2 ooo = +31 13'35.4", 
henceforth 4A) and NGC 1333 IRAS 4B (a 2 ooo = 3h29ml2.06s, <5 20 oo = +31° 13' 10.8", henceforth 
4B). Note that a 2 ooo and <5 2 ooo will denote J2000 epoch Right Ascension and Declination coordi- 
nates throughout this paper, respectively (a and 5 Q will denote offsets in Right Ascension and 
Declination from a p articular reference point, respectively). The s ource 4A has been extensively 



studied in the past (jSandell et al. 



1991 



Di Francesco et al 



2001 



Girart et al.1 120061 1. We have 



carried out continuum polarization and spectroscopic observations using, respectively, SHARP at 
350 /im towards 4A/4B and the 300 — 400GHz heterodyne receiver at the Caltech Submillimeter 
Observatory (CSO) on 4A. It is these data and the subsequent analysis that are presented. 

Previous work has greatly improved our knowledge of both 4A and 4B. The interferometric 
detection of inverse P-Cygni profiles towards both of these sources provides strong evidence for 
infalling gas motions (Di Francesco et al. 2001, hereafter DF2001). The detections were made in 
H 2 CO 3 2 i — > 2n emission and have allowed the determination of infall speeds of 0.68 and 0.47km/s 
for 4A and 4B, respectively. DF2001 also provide simple mass estimates for the gas in each source 
and find 0.71 M and 0.23 Mq within corresponding radii of 9"(0.013 pc) and 6" (0.009 pc) from 
the peak flux of 4A and 4B, respectively. Finally, mass accretion rates of 1.1 x 10 -4 M.@/yr and 
3.7 x 10~ 5 M /yr were calculated for both 4A and 4B, respectively. 



Recent 877 /im continuum polarimetry done at the Submillimeter Array (SMA) has indicated 



lr The distance of 300 pc will be assumed throughout this work. 
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the presence of a well denned pinch in the magnetic field morphology around 4A (Girart et al. 
2006, hereafter G2006). From these measurements the authors estimate the mass-to-flux ratio to 
be ~ 1.7 times the critical value of collapse. Along with their computation of the ratio of the 
turbulent to magnetic energy, /?turb = ^turb/^A ~ ^ x 10~ 4 (<5$int/l°) 2 = 0.02 (where o" tur b is the 
turbulent line width, Va is the Alfven speed, and 59- mt is the intrinsic dispersion in the polarization 
vectors; Lai et al. 2002), the authors conclude that 4A is an example of a magnetically dominated 
collapsing cloud core. A cloud mass of 1.2 within a radial distance of 3" (0.004 pc) from the 
peak flux of 4A is also traced by their dust continuum measurements. Taken together, the DF2001 
& G2006 results present compelling evidence for the notion that the physics of 4A is consistent 
with standard magnetically-regulated star formation theory. 

One problem not addressed in these works is the variation of physical parameters as a function 
of spatial scales. Of particular interest here is the variation in the magnetic field morphology with 
spatial scales; models predict that regions undergoing collapse will drag in the field lines towards 
the central condensation producing an hourglass morpholo gy. Further out from the con densation, 



the field morphology should remain in its ambient state (IFielder Sz Mouschoviad Il993l ). The re- 
sults presented in G2006 illustrate an example of this hourglass morphology at a resolution of 
~ 1" (0.001 pc). However, due to this small spatial scale G2006 was unable to sample the field 
morphology at scales larger than ~ 10" (0.015 pc) where models predict the field to be uniform. 
Similarly for the spectroscopy work, the high resolution attained through interferometry by DF2001 
(« 2" or 0.003 pc) allowed them to identify infall signatures out to a distance of ~ 4" (0.006 pc) from 
the peak flux of 4A and 4B. It is not clear from their results whether or not infall motions are 
occurring further out from the peak positions. 

The aim of this study is to address the problem of spatial scale variation in magnetic fields and 
infall motions and to complement the work of G2001 and DF2001. This will be done by analyzing 
single-dish observations obtained at larger spatial scales and comparing these results with the 
aforementioned papers. In the following sections we describe both polarimetric and spectroscopic 
observations carried out at the CSO and discuss the implications of our findings. In Section [2] we 
discuss our observations. Section [3] will cover a general discussion and analysis of our results, and 
finally in Section [4] we state our conclusions. 



2. Observations 

The following two sub-sections describe the submillimeter dust continuum polarimetry and 
spectroscopic observations acquired by our group. In Section 12.11 we discuss 350 \im polarimetry 
data collected in September 2008 using SHARP, the SHARC-II polarimeter. Section ^. 2l will discuss 
the HCN J = 4 — ► 3 spectroscopic observations taken with the CSO 300 — 400 GHz heterodyne 
receiver in September 2000. 
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2.1. Polarimetry 



Dust continuum polarimetry was done with SHARP at 350 //m with a spatial resolution of 
10" (0.015 pc). SHARP is a fore-optics additi on to the SHA RC-II camera that enables this 
instrument to be used as a sensitive polarimeter (|Li et al.ll2008l ). Although both 4A and 4B are 
studied here, the telescope was pointed onto 4B so as to provide the best possible chance of detecting 
polarization on this fainter source. Our map was calibrated with approximately 3 hours worth of 
data obtained on W3(OH) during the same observing run. In order to account for any random or 
systematic uncertainties that may remain after applying our standard data reduction pipeline, a 
reduced-;^ 2 analysis was performed on the data. The measurement uncertainty was correspondingly 
inflated, on a pixel- by-pixel basis, such that the reduced-^ 2 = 1. Our results are shown in Figure 
1 and represent approximately 10.5 hrs of observing time. 

Several points are worth mentioning from Figure 1: (1) the extended magnetic field around 
4 A is clearly sampled in our map. The nature of the extended field appears to be uniform out 
to a distance of ~ 20" (0.03 pc) from the peak of this source. This is consistent with the 85 fim, 
effective 20" beamwidth, SCUPOL observations taken at the JCMT (jMatthews et alJbood ). Our 
data enables a rough upper-limit to be placed on the size of the magnetic pinch reported in G2006 
to one SHARP resolution element (~ 10" or 0.015 pc). (2) deviations in the field appear as one 
moves out beyond 20" (0.03 pc) from 4A towards 4B. In the vicinity of 4B, the field morphology 
is rotated by ~ 30° towards the horizontal with respect to the field orientation around 4A. (3) 
depolarization is observed towards the peak of 4B, as denoted by open circles on the map where 
p + 2<7 p < 1%. Here p is the degree of polarization and a p is the corresponding uncertainty in the 
value of p. Significant depolarization is not observed towards 4A. 

Taken together, these results suggest that 4A and 4B are magnetically distinct objects. While 
at first glance the polarimetry of 4A appears to be consistent with conventional ideas of magnetic 
support, 4B is a more complicated case. The observed depolariz ation on 4B may result from changes 



in th e dust grain properties or shapes due to grain-growth (jVrba et al.l Il993l ; iHildebrand et al 



1999 ), supersonic an d super- Alfvenic turbulence plus a lack of grain alignment above Ay ~ 3 mag 

magnetic field geometr y and the inclination angle with the line-of-sight 



Padoan et al 



2001 



Goncalves et al 



scale field structures (Rao et al 



2005j ^Ficgc fc Pudritd l2000l ). or finally because of beam smearing over small- 



1998). The la st point cou ld be tested in a straightforward manner 



with high-resolution polarimetry. In addition, IChoil (|200ll ) identify a Class I object (denoted 4BII 
in the paper) within close proximity to the younger 4B core. The existence of this more evolved 
object is linked to the "Cav2" cavity, located to the west o f the 4B co mplex, and may have gener- 
ated this feature through the action of an ancient outflow ( Choil I2OO1I ) . It is possible to speculate 
that turbulence driven by this ancient outflow activity from 4BII may have disrupted the embedded 
magnetic field and thus contributed to the observed depola rization. W e should note that no such 
Class I objects have been identified within the 4A complex (|Choill2005l ). 



We now look at the position angles of the vectors situated within 20" of the peak of 4A in 
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order to assess the orientation of the large-scale magnetic field around this source. We assume this 
large-scale field has a simple, uniform morphology with a polarization angle 9. The value of 9 can 
be calculated by computing the mean of the individual orientation angles 9i, where the subscript 
i denotes the individual vectors. Table [T] contains all the relevant information on all the vectors 
depicted in Figure 1. 



From the data presented in Table HJ a straight arithmetic average of the orientation angles 9i 
for the vectors situated within a 20" radius of the peak of 4A yields a mean value of 9 ~ 45.9°. 
The dispersion in the orientation angles for these vectors (50 o bs) is found to be « ±13.6°. We 
therefore adopt 9 ~ 45.9° as the orientation angle of the large-scale uniform magnetic field around 
4A. The G2006 result of ~ 61° for their field orientation is approximately a la deviation from our 
result, implying the two values are consistent with one another. Uncertainties in our methods could 
account for the ~ 15° discrepancy; we admittedly assume a very simple uniform model for the field 
orientation, ignoring any possible non-uniform large scale structure to the morphology. The main 
outflow from 4A extends ±2' with a position angle of about 45° for a line drawn from one tip of the 
outflow to the othe r, but this rotates to 19° for t he inner part of the outflow confined to a radius of 



approximately 40 // (|Blake et al.lll995l ; IChoil 120051 ) . The reason for the change in orientation of the 
outflow is unknown; cloud core rotation may have played a role. A discussion of the full implications 
of this finding on the role of rotational support in the formation of 4A is beyond the scope of this 
paper. One would expect to see a significant deviation in the fi eld morphology from small to large 



spatial scales if centrifugal forces were dominant in this system (jMachida et al.ll2006l ). Instead, the 
polarimetry results presented here plus those obtained at smaller spatial scales (G2006) indicate 
otherwise. Nevertheless, the presence of a binary system within 4A shows centrifugal forces could 
not be negligible in the formation of this system. All of this information is consistent with the idea 
that the magnetic and centrifugal forces were comparable in magnitude for this system during the 
onset of collapse (G2006). In addition, it should be noted that the large-scale uniform magnetic 
field implied by our results is aligned with the original (i.e., large-scale) orientation of the outflow 
from 4A. 

Finally, we wish to calculate the mean intrinsic dispersion angle (£#i n t) over 4A by comparing 
our observations with the adopted uniform magnetic field model with 9 ~ 45.9°. Now 59^ is given 
by S9f nt = 59^ hs — a^, where cjq is the uncertainty in the observed position angle. We calculate 
values of <5# bs ~ 13.6° and a$ ~ 7.7° and as such work out 59- mt to be 11.2°. We will employ this 
value in our general discussion in Section [3l Note that we do not attempt a similar analysis with 
our results over 4B. It is apparent from Figure 1 that we cannot fit our vectors over this source to 
a simple model for the field orientation (plus no polarization is detected over the peak flux of this 
source). As such, it is not possible to calculate a meaningful average field orientation 9 or intrinsic 
dispersion <56*i n t for 4B. 
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2.2. Spectroscopy 

Observations of the HCN J = 4 — > 3 (354.505 GHz) rotational transition from 4A were made 
with the CSO 300-400 GHz heterodyne receiver. The beam size at this frequency is approximately 
20" (0.029 pc) and thus samples a far larger region of space than was obtained with DF2001. 
Detections were obtained for a sequence of points lying approximately along the outflow axis of 4A, 
as well as for a single point that was displaced from the center in the perpendicular direction. In 
total seven different positions were looked at, the results are illustrated in Figure 2. 

What is immediately clear is the presence of a peak followed by a dip in the line emission 
centered on 4A and the two positions lying closest to it along the red-lobe of the outflow. We should 
note that the presence of the red-lobe outflow will distort our spectra, as the background that our 
source absorbs against is not flat but the outflow emission itself. Therefore an interpretation of 
our data becomes clearer once the outflow component of the spectra is fitted with a Gaussian and 
removed, as was done using the CLASS software package and is shown in Figure 3. The spectra 
in Figure 3 are characterised solely by the aforementioned emission peak and dip, with the dip 
being situated to the right of the peak in each case. Each maxima and minima are observed at 
velocities of 6.8 km/s and 8.0 km/s, respectively. The full-width-half-maximum of each peak and 
dip profile is approximately ~ 1.5 km/s, and the point between the maxima and minima where the 
line temperature Tb equals ~ K corresponds to a velocity of approximately V ~ 7 km/s in each 
case. These data have exactly the same characteristics as the H2CO 3i2 — ► 2n spectra presented in 
Figure 4 of DF2001 and are thus characteristic of inverse P-Cygni profiles. With this interpretation 
we note that the spectral signatures seen in Figure 2 are characteristic of an infalling envelope of 
gas around 4A plus a outflow. These features are not observed towards the blue-lobe of the outflow. 
This is to be expected due to the position of the blue-lobe outflow, which is situated in between 
the infalling envelope and the observer. Therefore, this component of the outflow does not provide 
a background against which the infalling material can absorb. We do not have spectroscopy data 
on 4B at this time, and as such we cannot comment on the nature of gas motions around this core. 

One also notices the disappearance of the inverse P-Cygni profile at a offset position of a Q = 
27", 6 = 40" along the red-lobe outflow. We can therefore state that the peak of 4A is surrounded 
by an infalling envelope of radius r c on the plane of the sky, where 0.05 pc < r c < 0.07 pc 
. The lower bound of r c is provided by the spectra obtained at an offset position of a Q = 15", 
5 = 30"; the furthest position away from the peak of 4 A at which the infall signature is still 
clearly apparent. The range of r r specified here is consiste nt with the value of ~ 0.1 pc predicted by 



models of self-gravitating cores (jMcKee &: Ostrikeii 120071 ). Note that this estimate relies upon the 



assumption that nothing obstructs our view of the infall signature at a radial distance ~ 0.07 pc 
from the peak flux of 4A. Adopting r c ~ 0.06 pc, we immediately note that the size of the infalling 
envelope is approximately ~ 4 times larger than the upper limit placed on the size of the magnetic 
pinch observed around 4A (see Sect ion 12.11). This resu lt is consistent with theoretical work on the 



collapse of magnetized cloud cores (jGalli &: Shulll993l ). 
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3. Discussion 
3.1. Kinematics 

To further our investigation of the 4A system we fit our three inve rse P-Cygni profile detections 



to an enhanced version of the "two-layer" model originally devised bv lMvers et al.l (|1996l ) and later 
used by DF2001. In its most general form, our treatment of this problem is to envisage two parallel 
layers of material moving towards an opaque central condensation. These two layers, denoted "front 
slab" and "rear slab" , represent the infalling envelope. The opaque source is taken to fill a fraction 
T of the telescope beam. This setup is the model of DF2001, where the presence of outflows was not 
accounted for. For our model, two more layers of material are included in this system to represent 
the blue- and red-lobe outflows. This scenario is illustrated in Figure 4. 

Quantifying this model into an expression of the line brightness temperature Tq (V) at a certain 
velocity V we obtain equation (1) below. Note the subscripts l B\ l R', '/', V, 'c', and '&<?' represent 
the blue-lobe outflow, red-lobe outflow, front slab, rear slab, central source, and cosmic background 
parameters, respectively. Therefore the analytical model for the scenario depicted in Figure 4 is: 



Tb(V) 



T) ■ J {T R ) ■ (1 - e- T *) ■ e- T ° + { 
+J(T / )-(l-e~ T /)-e 
-T • J(T C ) • (1 - e- T f- TB ) - 



r + (l _ T 

TB + J{T B ).(\ 
{l-T).J{T bg 



J {%)■(! -e 

■e~ TB ) 
(l-e 



Tr ) • e ~ T f~ TB 



(1) 



where J (T) = Tq/ (e T °/ T — l) is the Planck temperature as a function of the blackbody temperature 
T, and To = hv/k where h is the Planck constant, k is the Boltzmann constant, and v is the 
frequency. The optical depth for model component i is den oted as 73. Finally we note that tq = 
r r + Tf + tb and t* = tr + tq . Following iMvers et al.l (|1996l ) we model the different optical depths 
with Gaussian profiles: 



n = m ■ exp 



tqj ■ exp 



(v-Vi- v LSR y 

2a} 



(V + Vj - V LS rY 
2a] 



(2) 



(3) 



where i = f, R and j = r, B, r$i and Toj are the peak optical depths, a\ and aj are the respective 
velocity dispersions, and the velocity for the local standard of rest is taken to be Vlsr ~ 6.96 km/s 
(DF2001). Note that Vt and V r are the infall velocities for the front and rear slabs, respectively. 

We can simplify (1) by making note of two important properties of our data. First, our large 
beam size implies T ~ (DF2001 get a value of T ~ 0.3 with a 2" beam). Second, we neglect 
the contribution of the blue-lobe, as only the red-lobe component of the outflow provides the 
background radiation against which the infalling material can absorb. It is clear from Figure 2 that 
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the spectra centered on (a Q = 0", 5 Q = 0") has a significant blue-lobe outflow component while the 
spectra centered on (a D = 7", 5 Q = 20") and (a Q = 15", 5 Q = 30") are largely dominated by the 
red-lobe outflow. As such, our outflow approximation will be especially coarse in the case of the 
spectra centered on the peak flux of 4A. By setting T = and tb = in equation (1) we get: 



T B {V,T B = 0) 



(4) 



J (T R ) • (1 - e"**) • e- T f-* + J (T r ) • (1 - e"^) 
+J(T / )-(l-e-/)- J(T bg )-(l-e-^), 

where % = tj + r r + tr. A computer program was written to minimize the reduced-x 2 function 
generated by comparing equation (4) with a multi-Gaussian fit to the three aforementioned spectra 
in Figure 2 that exhibit inverse P-Cygni profiles. T he minimization w as carried out through the use 
of Powell's Method for multidimensional functions ([Press et al.ll2002l ). The resulting fit of equation 
(4) to our spectra at offset positions (a = 0", 5 Q = 0"), (a Q = 7", 5 Q = 20"), and (a a = 15", 
5 = 30") are shown in Figure 5. 

The plots shown in Figure 5 demonstrate a reasonably good agreement between the actual data 
and the model. We should take stock at this point to stress that this model provides the simplest 
possible explanation for a contracting system with outflow. As such, it will provide us with only an 
approximate picture of the physical properties at play in 4A and its surroundings. This is especially 
true with regards to the optical depths and temperatures of the absorbing/emitting gas (DF2001). 
It should also be mentioned that the resolution of our data prevents us from distinguishing between 
small-scale structures, such as the binary system at the core of the 4A complex (dubbed 4A1 and 
4A2, G2006). Although we dispense with the optically thick c entral cond ensation in equation (4), 
each member of the binary system drives a unique outflow ( Choi 2005l j 2 1. By far the dominant 
outflow originates from 4A2, which po ssesses a length ~ 11 times l onger and is more luminous 
than the outflow associated with 4A1 (jChoil 120051 : iBlake et al.lll995l ). Our spectra likely sample 
an average of the environment of 4A, with detections of both outflows being integrated for data 
obtained on the peak flux of 4A. However, it is unlikely the data obtained away from the peak will 
be affected by the outflow of 4A1, since this outflow is limited in extent and has a position angle 
shifted by ~ 20° to the position angle of the much larger 4A2 outflow. Despite these shortcomings, 
our model is still useful for the comparison of the infall velocity between different spectra, since the 
ratio Vf/<7f (or V r /a r ) is a key parameter for any model of a contracting system (DF2001; Leung 
& Brown 1977). 

The model parameters resulting from the simulations shown in Figure 5 are listed in Table 
[2j One notices immediately that Vf > V r for each scan, where the values of Vf tend to be close 
to a speed of ~ 1 km/s while V r tends to have values closer to ~ 0.3 km/s. A likely explanation 
for this is that the value of Vlsr set in our model is somewhat inaccurate, and thus introduces 
a systematic error in our calculations. We note DF2001 leave Vlsr as a free parameter in their 
version of the "two-layered" model. Our attempt to treat Vlsr as a free parameter in equation 



2 Note that all referrals to the "4A outflow" made elsewhere in this paper pertain to the one driven by 4A2. 
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(4) failed to produce satisfactory results. This may be due to the highly non-linear nature of our 
model, which makes an iterative fit of equation (4) to a particular data set very sensitive to the 
initial parameter values. The simplest means to proceed was by fixing Vlsr to the DF2001 value. 
The differences between equation (4) and the DF2001 model may account for the different Vlsr 
value in our case. To correct for this we take an arithmetic average of the Vf and V r values listed in 
Table El which result in mean infall velocities of 0.63 km/s, 0.61 km/s, and 0.67 km/s for the scans 
at offset positions (a Q = 0", S = 0") , (a Q = 7", 5 Q = 20"), and (a Q = 15", 5 a = 30"), respectively. 
Taking all three scans together we estimate a mean infall velocity of ~ 0.64 km/s for 4A, a value 
that is very si milar to the 0.68 km /s velocity calculated by DF2001. If we assume a gas temperature 



of T g « 30 K (jBlake et al.l 1 19951 ) we can calculate the isothermal sound speed V ims = y/kT g /fimu 
to be ~ 0.33 km/s, where k is Boltzmann's constant, \x is the mean molecular weight (/z = 2.22), 
and tbh is the mass of hydrogen. Thus the observed infall is also supersonic. 

Assuming inside-out collapse for 4A, we note that the radius r c of the expansion wave for the 
infalling envelope is given by r c = V rms x t, where Vr ms is defined above and t is the time since the 
onset of collapse. Taking r c ~ 0.06 pc (see Section \2.2h . we find t ~ 2 x 10 5 yr. A large degree of 
uncertainty exists with regards to the age of the outflow, but estimates range from 2000 — 20000 yr 



(jBlake et al.lll995l ). Therefore the age of the infalling envelop in 4A is roughly 10 — 100 times the 



age of the associated outflow. 

Outflow velocities along the line-of-sight of —2.86 km/s, 10.21 km/s, and 9.17 km/s are esti- 
mated from the fits of equation (4) to the scans at offset positions (a = 0", 5 a = 0") , (a D = 7", 
5 Q = 20"), and (a Q = 15", 5 a = 30"), respectively. The latter two values indicate the red-lobe 
outflow has an approximate velocity of ~ 10 km/s. The negative value for the velocity at position 
(a Q = 0", 5 a = 0") is to due to the fact that equation (4) treats for a red-lobe outflow only, while 
the spectrum at this position shows a significant blue-lobe outflow as well. The velocity values for 
both red- lo be outflows are approximately consistent with SiO J = 1 — > observations presented in 



Choil (|2005l ) and illustrated in Figure 5 of that paper. The outflow observed at position (a Q = 0", 
5 Q = 0") is a composite of both the red- and blue-lobes and hence our model does not provide a 
reliable outflow velocity for this case. 

As a final note we wish to briefly discuss our handling of the outflow velocity Vr employed in 
our model for the scan at position (a a = 0", 5 a = 0"). Because of the simplifications introduced 
into equation (4), where only the red-lobe outflow is included analytically, it is necessary to set 
this parameter to a negative value in order to fit the composite red- and blue-lobe outflow that are 
present in this spectrum. Although this points out a shortcoming of our model, the important point 
here is to maintain emitting gas in the background of the infalling envelope; a red-lobe outflow. This 
component of the model provides emission from a backdrop of material against which the infalling 
gas can absorb and thus produce the dip present in a inverse P-Cygni profile. This is the reason 
why Vr is negative in equation (4) at position (a Q = 0", 5 = 0"). Despite this obvious defect, 
the model used here is the simplest mathematical construct that still conveys physical meaning for 
the data at hand. Future work may wish to employ more sophisticated Monte Carlo techniques to 
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generate models that more accurately describe this source. 



3.2. Support Mechanism 

We now proceed to calculate a mass estimate for 4 A and 4B. This can be achieved through a 
careful treatment of the thermal dust continuum data from SHARP. The mass of a cloud can be 
estimated via its submillimeter thermal emission by: 



(5) 



where R is the gas to dust ratio, g is the density of the dust material, d is the distance to the 
cloud, F v is the measured flux density, B v (T<j) is t he Planck function with dust temperature Td, 
a is the grain radius, and e u is the dust emissivity (|Hildebrandlll983l ). Here we assume values of 
R = 100, g = 3 g/cm 3 , d = 300 pc, a = 10" 5 cm, e v = 2.26 x 10~ 4 , and T d = 50 K, where the dust 
temperature was chosen to be the same v a lue as that employed in G2006. The values of g, a, and 
e v were taken or inferred from iHildebrandl (|1983l ). The total amount of 350 pm. flux detected within 
a radial distance of 20"(0.03 pc) from 4A is ~ 179 Jy, while a total of ~ 76 Jy was detected within 
10"(0.015 pc) of 4B. The distance around 4A was selected to encompass the polarimetry vectors 
used in the calculation of 56[ n t in Section 12.11 With this information, the mass of 4A and 4B is 
estimated to be « 1.9 M and « 0.8 M , respectively. These values are comparable to the mass 
estimates of DF2001 and G2006 that are stated in Section [T] of this work. As a final point we note 
that the total flux detected in Figure 1 is ~ 511 Jy. From this flux value a total detected mass of 
~ 5.4 Mq is calculated for the area of NGC 1333 IRAS 4 mapped by our continuum observations. 



We can now apply a modified form of the Chandrasekhar & Fermi (CF) technique (jChandrasekhar fc Fermi 



19531 ) to obtain an estimate of the magnetic field strength around 4A. If we assume the dispersion 
in our magnetic field map (see Figure IB) is entirely due to Alfven waves and/or turbulence, then 
the plane-of-sky magnetic field strength will be given by: 



B 



pos 



Q ( fe ) (4vrp) 



,1/2 



(6) 



Here we take Q 

of turbulent molecular clouds (lOstriker et al. 



0.5, the sa me value used in G2006 and one that is indicated by simulations 

We have already calculated 59 Ui t ~ 11.2° in 



2001 



Section [2j The gas density p can be roughly estimated by dividing the mass of 4A calculated above 
by the volume of a 20" (0.03 pc) radius sphere. Doing this yields a density of p ~ 8 x 10 -19 g/cm 3 . 
The line-of-sight velocity dispersion 8v\ os is taken from observations of HCO + J = 4 — » 3 obtained 
at the same time (and using the same heterodyne receiver) as the HCN observations discussed 
earlier. These data were obtained while pointing on the peak flux of 4A and are shown in Figure 
6. The choice of analyzing this spectrum for the 8v\ os value is motivated by the fact that an ion 
will be better coupled to the magnetic field (and the dust) than a neutral species over the whole 
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turbulent energy density spectrum (|Li Houdeil2008l ). Thus the line width of this species will be 



more representative of the turbulence/Alfven waves that may be disturbing the field as opposed 
to a neutral molecular counterpart. The value of 5v\ os is found to be 1.67 km/s. Inserting this 
information into equation (6) yields a field strength of -B pos ~ 1-4 mG, a value that is roughly a 
factor of 3.6 lower than the G2006 result. This difference is not unexpected since the field strength 
should increase towards the center of the core, as this is the location of maximum compression of 
the field lines. One could thus expect the value of -B pos to diminish at larger spatial scales. 

We are now in a position to calculate the mass-to-flux ratio for our observations on 4A. Writing 
this quantity in terms of the critical value for collapse A we find: 



A 



(7) 



where c$ = 0.12 (|Tomisaka et al.lll988l ). <J> is the magnetic flux, and G is the gravitational con- 
stant. Using this relation we calculate A ~ 0.44 over a region of ~ 270 arcsec 2 centered on the 
peak of 4A. This value is a factor of ~ 2.3 lower than unity and also notably lower than the G2006 
result. A strongly sub-critical cloud is also inconsistent with our own observations of infalling gas 
motions (Section 12. 2[) . Uncertainties in our estimate of M may be partly responsible; submillime- 
ter continuum emission will fail to sample the hot protostellar mass where dust particles cannot 
exist. We should note that the use of the plane-of-sky component of the magnetic field will result 
in an underestimate of and thus lead to a overestimation of A. Recent work suggests the incli- 
nation angle $ of the embedded magnetic field in this cloud falls within the range 0° < i? < 55° 
(jGoncalves et al.ll2008l ). This suggests that the magnitude of B = B pos /cos$ could fall somewhere 
within 1.4 mG < B < 2.4 mG. Therefore we can speculate that this effect could result in an 
overestimation of A by up to a factor of ~ 1.7. 

A more likely explanation for the sub-critical value of A may come from the application of 
the CF technique itself. Crucial to this calculation is an accurate measure of the intrinsic disper- 
sion angle 56- m t of the magnetic field vectors. Polarimetry, however, does not distinguish between 
contributions along the line-of-sight; the result being that the angular dispersion will be reduced 
through the process of signal integratio n through the thickness of the cloud as well as across the 
area subtended by the telescope beam (jHildebrand et al.l 120091 ) . The value of 59 mi that is then 
calculated could be smaller than the true dispersion and thus may result in an overestimation of 
-Bp OS through equation (6). Despite the fact that the fa ctor Q in equation ( 6) is meant to account 



for the problem of signal integration through the cloud (jOstriker et al.ll200ll). in reality th is is only 
a first-order correction and the smoothing effect could be more severe (jHoude et al.l 12009 ) . Higher 
resolution polarimetry across the large spatial scales observed here in conjunction with a more in- 



depth analysis may be used in the future to resolve this issue (jHildebrand et al. 
20091 ). 



2009 



Houde et al 



We can now calculate the ratio of the turbulent to magnetic energy (/?turb) i n 4A. The value 
we determine for this quantity is /?turb ~ 8 x 10~ 4 {56- in t/l°) 2 = 0.05, suggesting magnetic forces 
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dominate the physics of this particular cloud. From our discussion above we note that the field 
strength calculated here could be overestimated by a factor of a few. This could result from an 
underestimate of the intrinsic dispersion in the polarimetry 59- irA . For example, a requirement for a 
critical value of A would be an intrinsic dispersion value of 50- m t ~ 26°. This would then correspond 
to a value of /3 tur b ~ 0.5, a value sufficiently close to unity to suggest an equipartition of the 
turbulent and magnetic energies in this system. These results are definitely not consistent with 
is the idea of a turbulent-dominated support system for this cloud. Since models of magnetically- 
supported star formation predict A w 1 for cloud cores undergoing collapse, we conclude that our 
results are likely consistent with the idea that the turbulent and magnetic energies are of the same 
order of magnitude in this system. 



4. Conclusions 



We summarize our results as follows: 



• Dust continuum polarimetry done at 350 /iin with SHARP demonstrates a uniform magnetic 
field morphology around 4A at a resolution scale of 10" (0.015 pc). We therefore adopt the 
size of one resolution element to be an approximate upper limit on the size of the magnetic 
pinch region reported in G2006. This is in agreement with magnetic cloud collapse theory 
and correlates well with the findings of G2006. In addition, the large-scale uniform magnetic 
field appears to be aligned with the original (large-scale) outflow direction for 4A. Why the 
large-scale and small-scale directions of the 4A outflow are different is not known, but our 
results and those of G2006 show that the average orientation of the magnetic field has not 
changed direction from large-scales to small. The polarimetry obtained on 4B appears to 
indicate depolarization towards the peak flux region. An explanation for these observations 
on 4B remains an open point. Mass estimates for both 4A and 4B have been made revealing 
1.9 M within 20"(0.03 pc) and 0.8 M within 10" (0.015 pc) of the peak flux for both cores, 
respectively. The total mass traced by our continuum observations is 5.4 M for the portion 
of the NGC 1333 IRAS 4 cloud complex surveyed. 

• Spectroscopy done at the CSO with HCN J = 4 — > 3 line emission has revealed inverse P- 
Cygni profiles at offset positions (a Q = 0", 5 Q = 0"), (a = 7", S Q = 20"), and (a = 15", 
5 = 3 0") from the l ocatio n of peak flux for 4A. Fitting these data with an enhanced version 



of the I Myers et al.1 (|1996l ) "two-layer" model, we estimate a mean infall speed of 0.64 km/s 
for this cloud core. These findings are in good agreement with the results of DF2001. The 
radial size of the infalling envelope r c is estimated to range between 0.05 pc < r c < 0.07 pc 
and is thus approximately ~ 4 times larger than the size of the magnetic pinc h. This is 



consistent with theoretical work on magnetized cloud collapse (jGalli Shulll993l ). The age 
of the infalling envelope is found to be approximately ~ 2 x 10 5 yr; a figure that is roughly 
10 — 100 times larger than the age of the associated outflow. 
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• The value of the plane-of-sky component of the magnetic field strength has been estimated 
to be ~ 1.4 mG around 4A. This yields a normalized mass to flux ratio of A ~ 0.44. This 
value of A is inconsistent with our observations of infalling gas motions and does not agree 
with the previous results of G2006. We speculate here that <5#i n t may be underestimated due 
to smoothing effects in the angular dispersion along the line-of-sight that is characteristic of 
polarimetry. If we inflate our value of 59 ia t by a factor of 1/A , we then find that /?turb ~ 0.5. 
This suggests an equipartition of the magnetic and turbulent energies in this cloud. 
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Fig. 1. — SHARP polarimetry (a) and deduced magnetic field orientation (b) over 4A and 4B. 
Both images are centered on 4B («2000 = 3/i29ml2.06s, ^ooo = +31°13'10.8") with 4A lying 
towards the northwest corner of the map. The horizontal and vertical axes show offsets in Right 
Ascension and Declination, respectively. Contour levels are 0.1, 0.2, . . . 0.9 times the peak flux 
value (29.3 Jy /beam). Image b) also shows the G2006 magnetic field map for comparison, where 4A 
is resolved into its components 4A1 and 4A2. Arrows indicate the orientation of the outflow. Note 
that all the vectors presented in image a) are p > 3a p and circles denote regions where p+2a p < 1%. 
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NGC 1333 IRAS 4A - HCN(4-3) 
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Fig. 2.— HCN J = 4 -> 3 line emission (354.505 GHz) centered on 4A (a 2 ooo = 3/i29ml0.42s, 
£2000 = +31°13 / 35.4 // ). Observations were made out to offset positions of (a = —15", S = —30") 
on the blue lobe of the outflow and (a = 27", S = 40") on the red lobe. One observation at an 
offset position of (a Q = —30", 5 Q = 15") was made perpendicular to the outflow axis. 
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NGC 1333 IRAS 4A - HCN(4-3) 
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Fig. 3.— HCN J = 4 -> 3 line emission at offset positions (a a = 0", <5 = 0"), (a a = 7", S = 20"), 
and (a Q = 15", 5 = 30") from the peak flux of 4A. Outflow components have been removed 
to reveal characteristic inverse P-Cygni profiles at each of the three locations. The outflows were 
removed with the use of the CLASS software package. Vertical lines represent velocities at 6.8km/s 
and 8.0 km/s. 
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Fig. 4. — Illustr ation of the enhanced "four-layered" version of the two-layered model of 



(jMyers et al.lll996l ). This depiction labels each component of the model with a Planck temper- 
ature J(Tj), optical depth n, and a velocity The subscript i denotes l B\ l R\ '/', 'r', and 
'c' representing the blue outflow, red outflow, front slab, rear slab, and central source parameters 
respectively. Arrows indicate velocity directions along the observers line-of-sight. We also allow for 
cosmic background radiation with a Planck temperature J (Tft s ) = 0.49 K. 
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Fig. 5. — HCN J = 4^3 line emission observed at offset positions (a = 0", S Q = 0") (a), (a Q = 7", 
5 = 20") (b), and (a Q = 15", 5 Q = 30") (c) from the peak continuum flux of 4A. The dashed curves 
show the fit of our model, equation (4) in the text. 
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Fig. 6. — HCO + J = 4 — > 3 line emission at an offset position of (a = 0", <5 = 0") from the 
peak continuum flux of 4A. These data were used for the evaluation of Sv\ os in our calculation of 
Bpos- We have assumed the ion will be more closely associated with the magnetic field through 
the Lorentz force, and hence the line width observed here should be more representative of the 
turbulance/MHD waves that may be distrubing the embedded field. The line width is calculated 
to be 5v\ os ps 1.67 km/s. 
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Table 1. Measured Polarization for NGC 1333 IRAS 4 C . 
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a Note 6i angles describe the orientation of the 
deduced magnetic field. Angles are measured rel- 
ative to north and increasing eastward 

b Offset positions with respect to the peak posi- 
tion of NGC 1333 IRAS 4B 

c Data given below the solid line describe vectors 
associated with 4A 



Table 2. Model Parameters for HCN J = 4 -> 3 Inverse P Cygni Profiles 13 . 
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a Offset positions with respect to the peak position of NGC 1333 IRAS 4A 
b Variable definitions are given in Section 3.1 



